Tracer density discontinuities in turbulent flows: simple model and scaling laws 
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Mixing in fully developed incompressible turbulent flows is known to lead to a cascade of discon- 
tinuity fronts of passive scalar fields. A one-dimensional (ID) variant of Baker's map is developed, 
capturing the main mechanism responsible for the emergence of these discontinuities. For this ID 
model, expressions for the height-distribution function of the discontinuity fronts and structure 
function scaling exponents C,p are derived [for Kolmogorov turbulence, — |log3(p+ 1)]. These 
analytic findings are in a good agreement with both our ID simulations, and the results of earlier 
numerical and experimental studies. 
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Explaining the origin of intermittency in turbulent me- 
dia is a long-standing challenge for statistical physicists; 
a particular attention has been paid to the anomalous 
behavior of the structure functions scaling exponents C,p. 
While the first studies in this field date back to almost 
five decades 1], the theory is still far from being com- 
plete. Even in the simplest case of passive scalar tur- 
bulence, despite of extensive studies (c.f. reviews [1,0]), 
the theoretical understanding of the anomalous scaling is 
rather sketchy: the exponents are either obtained experi- 
mentally or numerically (e.g. [1, i, i, 0, S, i, El ) ; 

the an- 
alytic results are limited to specific velocity spectra (not 
applicable to the Kolmo gor ov case) 0, hj. There are 
also hierarchical models 1131 (derived from the velocity 
field intermittency model [l^), which can fit the experi- 
mental data relatively well (although being inconsistent 
with the Obukhov-Corrsin result for ^2 but rely on 

a couple of phenomenological hypothesis; therefore, their 
usefulness in understanding the origins of intermittency 
is limited. 

It is known that if a passive scalar evolves in fully devel- 
oped turbulent flows, the scalar fleld becomes everywhere 
discontinuous: the discontinuity fronts of fractal struc- 
ture will emerge These fronts are the very reason 
for the anomalous scaling of structure functions. How- 
ever, little is known about the formation and statistics of 
them. 

In the first part of the Letter, we outline qualitatively 
the mechanism of the formation of passive scalar discon- 
tinuity fronts, and construct a ID model incorporating all 
the essential features of that mechanism. In the second 
part, we analyze the properties of our model theoretically, 
applying a non-rigorous scaling analysis, present the sim- 
ulation results, and compare them (as well as some earlier 
experimental and numerical results) with the theoretical 
scaling laws. 

/. Formation of discontinuities. In what follows, we as- 
sume that the tracer density 9{r, t) evolution is described 
by a simple diffusion equation. 



dtO + vVe 



fir,t), 



(but not zero). Finally, it is assumed that the velocity 
field obeys a power-law Kraichnan statistics. 
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where dij{r) is the non-constant part of Dij{r), Di — 
a constant, and f — the smoothness exponent, c.f. 
and references therein. So, we neglect the intermittency 
of the underlying velocity field. It should be noted that 
according to the experimental and numerical evidence, 
the passive scalar intermittency is stronger than the ve- 
locity field intermittency (e.g. characterized by greater 
anomality 2^2 — C4); c.f. [la]. Thus, one can expect that 
the former dominates over the latter, and the behavior of 
tracers in real (intermittent) velocity fields is very similar 
to that of in idealized Gaussian fields (for a numerical ev- 
idence, c.f. [IB])- Finally, note that owing to the robust- 
ness of our model, the assumption of delta-correlation in 
time will not be actually used; Eq. ([2]) is adopted for a 
starting point merely to simplify comparisons with other 
studies. 

Mixing effect of turbulent flows is most intuitively 
characterized by the growth of the distance between two 
tracer particles r: 



^ (Inr) (X ^ ([Inr]^) oc r« ^, 



(3) 



c.f. [13] • So, with a proper time unit, the distance dou- 
bling time is estimated as r w r^"^; the Kolmogorov 
scaling r ~ r^^^ is matched with ^ = 4/3. For the sake 
of simplicity, we consider two-dimensional (2D) geometry 
(generalization to the 3D geometry is straightforward). 

The formation of tracer discontinuities can be quali- 
tatively explained as follows. First, we decompose the 
velocity field into components of different characteristic 
space-scale of size a, Va{r,t) — f^^^^^^.^^v{k,t)e^''^dk, 

where v{k,t) is the Fourier component of v{r,t). The 
characteristic time-scale is defined as the time needed 
for the field Va{r, t) to transport a tracer particle to a dis- 



(1) tance of the order of a; according to Eq. ([3]), Tq 
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where v{r,t) is a turbulent velocity field, f{r,t) is a forc- 
ing, and the seed diffusivity k is assumed to be very small 



Suppose that initially {t — 0), there is a constant tracer 
gradient: 0{r,Q) = x. In Fig. 1, the regions 9 < 0, 
< < 1, and 1 < 9 are marked by black, gray, and 
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FIG. 1: Simplified scheme of the formation of tracer discon- 
tinuities. Characteristic time of eddies of size a scales as 
r ~ a^~^. Due to the combined efi^ect of large and small 
eddies, low- and high-density regions (black and white, re- 
spectively) are brought into contact within a finite time. 



white, respectively. Minimal distance between the iso- 
Unes = and 9 = \ will start decreasing (they are 
turbulently stretched, and the surface are between them 
is conserved). Initially, this approaching is dominated 
by the largest eddies fitting between these lines, i.e. the 
eddies of approximately unit diameter, characterized by 
Ti « 1. Then, after a unit time, the distance between 
some segments of the isolines will be decreased approx- 
imately by a factor of two. From now on, the distance 
decreasing rate will be dominated by twice smaller ed- 
dies, and the characteristic time-scale is reduced by a 
factor of 22-«. The process will continue ad infinitum, 
leading to the contact of segments within finite time (the 
characteristic time scales form a geometric progression). 

We aim to construct a model, which mimics the evolu- 
tion of the tracer density profile along the x-axes (any ID 
cross-section). To begin with, we consider only the effect 
of an "a-flow" ' Va[r, t) (this corresponds to observing the 
initial tracer field evolution with a spatial resolution a: 
smaller vortices are not resolved, larger ones are slower 
and require a longer observation period). In incompress- 
ible velocity fields, exponential growth of scalar density 
gradients is due to exponential stretching of fluid ele- 
ments, caused by stretching-folding motion of the fluid 
(c.f. jlMl)- Such a stretching- folding motion is provided 
by a simple shear flow, as depicted in Fig. 2. Now, con- 
sider the tracer density profile along a;-axes in Fig 2: ini- 
tially monotonous curve 9{x) is replaced by a new profile 
9'{x) with a "kink" . The kink emerges because a descend- 
ing segment is substituted by the sequence of descending, 
ascending, and descending segments. In the idealized ver- 
sion, all these curve segments are mirror images of each 
other and the result of a 3-fold "compression" of the ini- 
tial curve segment along the cc-axes. Such a mapping 
Ma.c ■ 9{x) — > 9'{x) is represented in Fig. 3(A); here, a 
denotes the size of the vortex and c — its middlepoint. 

Our model for the tracer turbulence is iterative ap- 
plication of the mapping M.at,ct to some initial profile 
9tj{x) with random values of the parameters a„ and c„: 
9t+i{x) = Mat,ct [St{x)]] note that t plays the role of (dis- 
crete) time. In order to match statistically homogeneous 
turbulence, the probability distribution function (PDF) 
of this mapping over the parameter c needs to be homo- 
geneous, and PDF over a needs to match the stretching 
statistics of the velocity field Q. Therefore, the wait- 




FIG. 2: The effect of a shear fiow on the tracer density profile 
along the a;-axes. Suppose that initially, the tracer density 
gradient is constant (left image). Initially straight isolines 
evolve into s-like curves (black and white curves in right im- 
age; C and D are their touching points with the a;-axes). The 
segment AB of the initial profile 9{x) evolves into a "kink" 
of the final profile 6'{x), consisting of descending, ascending, 
and descending segments A'C, CD, and DB' , respectively. 




FIG. 3: Mapping A4a,c modeling the efi^ect of a single vortex 
of size a on the tracer profile 6{x). The vortex may be far 
from the vessel boundary [Case (A)], or it may be close to the 
boundary, at which a fixed tracer density value 9\x=o = 1 is 
kept [Case (B)]. 



ing time Ta between two subsequent mappings of size a 
at X = c (i.e. satisfying the conditions at G [a, 2a] and 
c e [ct — at/2, ct + at/2]) needs to scale as Ta ~ a^^^. 

We need also to address the issue of the boundary con- 
ditions. Initial conditions in the form 9{t = 0, r) are mod- 
elled by the initial profile 9o (x) . The situation when there 
is no tracer flux through the container boundaries can be 
modelled by periodic boundary condition 9t{x + 1) = 
9{x). The effect of forcing in Eq. ^ is mimicked by ad- 
ditional iterations 9t+i{x) = 9t(x) + f{t,x). If there is 
a boundary with a fixed tracer density 9\x — — 1, we 
need to incorporate a mechanism of tracer influx at that 
density. For real turbulent flows, that influx is provided 
by these vortices, which are in contact with the wall, and 
is proportional to the vortex size. A convenient way to 
match such a process is represented in Fig. 3(B). If the 
outer edge of the mapping A4at,ct happens to be beyond 
the container boundary x — (i.e. Ct < at/2), the pro- 
file 9t{x) is extended to the region a; < with the value, 
kept at the boundary. Then, the mapping can be ap- 
plied in the same way as described before [with a single 
modification: "compression" factor is increased from 3 
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FIG. 4: As a result of the mapping M2a,B, the old value of 
the mean density difference Aa(-B) defines the possible range 
of new values at trice smaller scale 6 = a/3: the smallest value 
is Ai,(F) = 0, and the largest one /^i,[E) = A6(G) = Aa(-B). 



to 3(| — Ct/at), so that the entire "vortex" will fit inside 
the region x > 0]. 

Second, we need to discuss the effect of seed diffusivity. 
For any non-zero diffusivity k, the diffusion smoothes the 
tracer density fluctuations at a microscale 5, for which 
the effective Peclet' number k, 1. From the equahty of 
diffusion and mixing times, Tg k, S^~^ « S'^ /k, we obtain 
5 « K^/^. fn order to take into account such a smoothing, 
the mapping A4at,ct is modified so that apart from the 
effect depicted in Fig. 3, it includes also averaging over a 
sliding window of width S. 

For numerical simulations, 6 serves as a natural dis- 
cretization step: tracer density profile is stored as an 
array Oi^t = OtiiS), i = 1,2,... TV with N = 
Then, the iteration formula is 6*^^4+1 = \ '}2,\k-j\<i^k,t, 
where j = ct — 3(« — Ct), if « — Ct £ [— |at, ^Ot], and 
j = ct + 3{i - Ct ± iflt), if i - Ct ± |at € [-^h, |at]- 
Here, at = at/ 5 and ct — Ct/5 are the discretized map- 
ping parameters. 

//. Scaling analysis of the model. Our scaling analy- 
sis is based on the probability density function (PDF) 
/a(A,) of the difference A„(x) = \ea[x + f ) - Oaix - f )| 
between the mean values of the tracer densities for neigh- 
boring segments of length a. Here, the local average is 
defined as Oa{x) = i j^^S 0{y)dy. 

To begin with, we consider, what will happen, if seg- 
ments AB and BC, characterized by mean densities 6ab 
and 0BC, are transformed by a mapping, see Fig. 4. The 
segment AB is transformed into three trice smaller seg- 
ments, two of which are marked as DE and GH; the 
mean densities for these segments are equal to that of the 
segment AB. The same applies to the segments BC, EE, 
and EG. So, the density drop between the segments DE 
and EE equals to that of between the segments AB and 
BC, i.e. to Aq(_B). However, there is no density drop be- 
tween the segments EE and EG. Apparently, the density 
drop A^/3(x) takes all the intermediate values between 
Aa{B) and 0, as x moves from E to E; the dependence 
on X is approximately linear (at least in the neighbor- 
hood of E). Consequently, as a result of the mapping, 
the probability fa{Aa)dAa, associated with the density 
drop Aa, contributes to the PDF fa/3 evenly over the 
range of values A^/s G [0, A^]. This mechanism relates 
all the values of fb{Ai,) with 6 = | to the values of fa{Aa) 
[because mappings are continuously being applied to the 




FIG. 5: Theoretical curve family of the cumulative density 
difference probabilities P — J^,-^^ fa{A')dA' is calculated 

according to Eq. (O for ^ = | and a = 5, 45, ... , 4^5 (solid 
curves). The grey dots indicate the corresponding data-series 
of our simulations. Insert: the theoretical (^p-curve (solid line; 
^ = |) is compared with various experiments and Navier- 
Stokes numerical simulations; the dashed line corresponds to 
the Kraichnan formula [l9l]. The data-series are from the 
following references: (o) — ll, (&) — [1 (c) — [ll, (rf) — 0], 

(e)-i, (/)-i, ig)-m 
profile 0{x)]; mathematically, 

/,(A)=^'/a(A')^ (4) 

(assuming that the maximal value of A is 1). 

It should be emphasized that Eq. ^ is obtained, us- 
ing two implicit assumptions, (i) We do not consider the 
effect of those mappings, the size of which is either signif- 
icantly larger or smaller than a and b. It can be argued 
that the effect of smaller size mappings is insignificant at 
our scale, because they preserve the average density 9h 
(this is true, if the mapping falls entirely into the seg- 
ment; if it falls at the edge, 9b will be changed, but the 
change remains relatively small). However, for very small 
values of ^ < ^Oj when small vortices are much more fre- 
quent than the large ones, this assumption will no longer 
be valid, (ii) We can neglect the effect of larger vortices. 
This is actually not true: larger-size mappings compress 
the profile without reducing the density drop A via the 
process depicted in Fig. 4. Such a process corresponds to 
a direct transfer /a (A) f^/siA), without the convolu- 
tion in Eq (jj]). So, in average, the profile will be com- 
pressed more than trice, before entering the convolution 
stage. Hence, the effect of larger vortices can be taken 
into account by using an effective, somewhat increased 
compression factor k = a/b > 3. 

Bearing in mind boundary conditions 0{0,t) = 1 
and 0{l,t) = 0, it is reasonable to assume that 
/i(A) = 1. Then, direct integration results in fa{Aa) = 
I ln(Aa)|"/n!, where n = — logj. a is an effective number 
of iterations. Now we can easily calculate the structure 
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FIG. 6: The dependence of the second order anomalous scal- 
ing exponent 2C,i — ^4 on the smoothness exponent ^. Solid 
line corresponds to the theoretical curve (2 — ^) logg |, dot- 
ted line — to the Kraichnan formula [l^, dashed line — to 
the perturbation theory ^2Cl| , open circles and stars — to the 
numerical results of 2D and 3D Lagrangian simulations [lH . 
and filled circles — to the simulations with our ID model. 

function scaling exponents C,p. Indeed, we expect that 
//a(A)A^'(iA oc a^p; the integral is easily taken, result- 
ing in (^p — logj.(p +1). Compa ring this expression with 
the classical result ^2 = 2 — ^ [i4| (which is valid both 
for tracer turbulence, and for our ID model), we obtain 
k = 3^/*^^"^-'. This equality allows us to rewrite the ex- 
pressions of /a and as 

/a(A) = |ln(A)l»n!-i, n = -(2 - ^ logg a, 
Cp = {2-e)log3(p+l). 

Now, let us recall that we expected fc > 3; this inequal- 



ity is not satisfied for ^ < 1. So, we can conclude that 
^0 = 1, i-e. for ^ < 1, the assumption (i) is not satisfied. 
Note that the result = 1 is directly applicable only to 
our ID model, when all the compression factors are equal 
to 3; in the case of real 2D or 3D turbulence, the effective 
compression factors may take different values and hence, 
the critical value may deviate from 1. 

We have implemented the above described model [with 
boundary conditions 9{0,t) = 1 and 0(1, t) = 0] numeri- 
cally for several values of £,. The array length was taken 
equal to TV = 10^ (recall that this parameter plays the 
role of the ratio of the tank width and diffusion scale). 
For each value of the observation time of the evolution 
of the density field dt,i was long enough to include at least 
10^ full decorrelations (i.e. at least 200 occurrences of the 
largest-sized mappings with a > Simulation results 
are presented in Fig. 5 and Fig. 6. The small mismatch 
between the curves and data points in left-hand-side of 
Fig. 5 can be explained by finite-size effects and some- 
what purer statistics of extreme events. The reason of 
the departure of the theoretical curve from the simula- 
tion data for ^ ^ 1 has been already discussed. 

In conclusion, our ID model and analytical results ex- 
plain, with a reasonable accuracy, the results of previ- 
ous experiments and simulations. Eq. ([5]) predicts that 
there is no saturation of the exponents Cp (for p — » oo). 
While some experiments have reported such a saturation 
(c.f. [§]), the others have not (c.f. [3]; the inconsistent 
results can be explained by pure statistics of extreme 
events (very large density differences) . The possibility to 
extend our approach to passive and active vectors (kine- 
matic dynamo and hydrodynamic turbulence problems) 
will be the scope of further studies. 

The support of Estonian Science Foundation grant 
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